# Figure 3
rm(list=ls())
gc()

# Load RCPs
## First column is temperature time series from deterministic EBM. initialised at 1C in 2020. Next 1000 columns are temperature time series from stochastic EBM. See "figure3.m" for details.
path <- "~/Documents/Projects/IAM_Stochastic/" # Directory where output from "figure3.m" is stored.
files <- list.files(path=path, pattern="output_")
d <- data.frame()
for (i in 1:length(files)){
	ds <- read.table(paste(path, files[i], sep=""), sep=",", header=FALSE)
	ds$rcp <- files[i]
	d <- rbind(d, ds)
}

# Ensemble size
nn = ncol(d) - 3
names(d) <- c("year","det",paste("X",seq(1:nn),sep=""), "rcp")
d$year <- d$year + 2019

# Length of time series
n=length(unique(d$year))

# Economic assumptions 
gdp = 80*10^12 # initial aggregate consumption (billion US dollars)
g = 0.019 # annual growth rate of undisturbed aggregate consumption
eta = c(0,1.45) # marginal utility of consumption
rho = c(0.04255,0.015) # pure rate of time preference (1.5% in DICE-2016)

s <- 1 # Index for assumptions about discount rate. When s=1, r = 4.255 + 0 * 1.9 = 4.255

# Demographic assumptions
pop = 7.5 # billion people
pop.asympt = 11.5 # Asymptotic population, in billion
param = 0.134/5 # DICE-2016 "Population 2050 parameter"
pop.ts <- rep(NA,n)
pop.ts[1] <- pop
for (i in 2:n) { pop.ts[i] <- pop.ts[i-1]*(pop.asympt/pop.ts[i-1])^param}
pop.ts <- matrix(rep(pop.ts,nn+1),ncol=nn+1,nrow=n)
pop.ts <- pop.ts*10^9

# Undisturbed output per capita
gdp.ts <- gdp*(1+g)^seq(0,n-1,1)
gdp.ts <- matrix(rep(gdp.ts,nn+1),ncol=nn+1,nrow=n)
gdp.ts <- (gdp.ts/pop.ts)/1000 # undisturbed consumption per capita expressed in thousands of dollars.

# Discount factors
df.2 <- 1/((1+rho[s])^seq(0,n-1,1))
df.2 <- matrix(rep(df.2,nn+1),ncol=nn+1,nrow=n)

x <- data.frame()
for (j in files) {
	t <- d[d$rcp==j,]
	t <- t[,!(names(t) %in% c("year","rcp"))]
	# Isoelastic utility
	x.Weitzman <- colSums((((gdp.ts*(1/(1+((t)/20.46)^2 + ((t)/6.081)^6)))^(1-eta[s]))/(1-eta[s])) * df.2 * pop.ts)
	x <- rbind(x, x.Weitzman)	
}
x <- t(x)
colnames(x) <- files

## Scaling factors
UofC <- ((gdp.ts[1,1])^(1-eta[s]))/(1-eta[s]) # Utility of undisturbed per capita consumption today (in thousands of dollars)
scalar <- UofC-(((gdp.ts[1,1]-0.001)^(1-eta[s]))/(1-eta[s])) # Utility value of the marginal dollar
UofC <- sum((((gdp.ts[,1])^(1-eta[s]))/(1-eta[s]))*pop.ts[,1]*df.2[,1]) # NPV of utility of undisturbed aggregate consumption

## Equilibrium Climate Sensitivity values for which temperature trajectories have been generated
ECS <- sapply(strsplit(colnames(x), "_"), "[", 2)
ECS <- sapply(strsplit(ECS, ".txt"), "[", 1)
ECS <- as.numeric(ECS)
ECS <- unique(ECS)

## Fit a log-Normal ECS distribution to IPCC most likely value and likely range.
fit <- data.frame(mean=seq(1.01,4,.01))
fit$sd <- (fit$mean - log(3))^.5 # Pairs each mean with a sd that produces a log-normal distribution with mode = 3.
fit$mode <- exp(fit$mean - fit$sd^2)
fit$dens.interval <- plnorm(4.5, meanlog = fit$mean, sdlog = fit$sd, log = FALSE) - plnorm(2, meanlog = fit$mean, sdlog = fit$sd, log = FALSE)
fitted.mean <- fit$mean[abs(.67-fit$dens.interval)==min(abs(.67-fit$dens.interval),na.rm=T) & !is.na(fit$dens.interval)]
fitted.sd <- fit$sd[abs(.67-fit$dens.interval)==min(abs(.67-fit$dens.interval),na.rm=T) & !is.na(fit$dens.interval)]
step <- max(ECS) - max(ECS[ECS<max(ECS)])
probs <- c()
for (i in seq_along(ECS)) {
	probs[i] <- plnorm(ECS[i]+step/2, meanlog = fitted.mean, sdlog = fitted.sd, log = FALSE) - plnorm(ECS[i]-step/2, meanlog = fitted.mean, sdlog = fitted.sd, log = FALSE)
}
# plot(ECS,probs)

## Truncates log-Normal distribution at ECS=10. Enable to calculate first row of SI table 1
# probs[ECS>10] <- 0
# probs <- probs/sum(probs)

## Fits a Roe/Baker distribution instead of a log-Normal. Enable to calculate third row of SI table 1
# meanF <- .65
# sigmaF <- .13
# T0 <- 1.2
# densT <- function(T) { (1/(sigmaF*sqrt(2*pi)))*(T0/T^2)*exp(-.5*((1-meanF-(T0/T))/sigmaF)^2) }
# step <- max(ECS) - max(ECS[ECS<max(ECS)])
# probs <- c()
# for (i in seq_along(ECS)) {
# 	probs[i] <- integrate(densT, lower = ECS[i]-step/2, upper = ECS[i]+step/2)$value
# }
# probs <- probs/sum(probs)
# plot(ECS,probs)


s <- 2 # Index for assumptions about discount rate. When s=2, r = 0.015 + 1.45 * 1.9 = 4.255

# Discount factors
df.1 <- 1/((1+rho[s])^seq(0,n-1,1))
df.1 <- matrix(rep(df.1,nn+1),ncol=nn+1,nrow=n)

x1 <- data.frame()
for (j in files) {
	t <- d[d$rcp==j,]
	t <- t[,!(names(t) %in% c("year","rcp"))]
	# Isoelastic utility
	x.Weitzman <- colSums((((gdp.ts*(1/(1+((t)/20.46)^2 + ((t)/6.081)^6)))^(1-eta[s]))/(1-eta[s])) * df.1 * pop.ts)
	x1 <- rbind(x1, x.Weitzman)	
}
x1 <- t(x1)
colnames(x1) <- files

## Scaling factors
UofC1 <- ((gdp.ts[1,1])^(1-eta[s]))/(1-eta[s]) # Utility of undisturbed per capita consumption today (in thousands of dollars)
scalar1 <- UofC1-(((gdp.ts[1,1]-0.001)^(1-eta[s]))/(1-eta[s])) # Utility value of the marginal dollar
UofC1 <- sum((((gdp.ts[,1])^(1-eta[s]))/(1-eta[s]))*pop.ts[,1]*df.1[,1]) # NPV of utility of undisturbed aggregate consumption


j = rev(files)[1]
(x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x),grep(j,colnames(x1))]))/scalar1 # Risk premium in present value dollars
((x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x),grep(j,colnames(x1))]))/scalar1)/gdp # Risk premium as a share of current global output

j = rev(files)[2]
(x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x),grep(j,colnames(x1))]))/scalar1 # Risk premium in present value dollars
((x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x),grep(j,colnames(x1))]))/scalar1)/gdp # Risk premium as a share of current global output

j = rev(files)[3]
(x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x),grep(j,colnames(x1))]))/scalar1 # Risk premium in present value dollars
((x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x),grep(j,colnames(x1))]))/scalar1)/gdp # Risk premium as a share of current global output

j = rev(files)[4]
(x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x),grep(j,colnames(x1))]))/scalar1 # Risk premium in present value dollars
((x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x),grep(j,colnames(x1))]))/scalar1)/gdp # Risk premium as a share of current global output


# # Unused plot. Displays the QQ-plots for the distributions of climate damages with and without aleatory uncertainty. Shows that they deviate primarily in the upper tail of the distribution, which is a consequence of the risk interaction effect discussed in the paper.
# ## IPCC AR5 colour scheme
# c <- data.frame(rcp=c("26","45","60","85"))
# c$r <- NA
# c$r[grep("26",c$rcp)] <- 0
# c$r[grep("45",c$rcp)] <- 121
# c$r[grep("60",c$rcp)] <- 255
# c$r[grep("85",c$rcp)] <- 255
# c$g <- NA
# c$g[grep("26",c$rcp)] <- 0
# c$g[grep("45",c$rcp)] <- 188
# c$g[grep("60",c$rcp)] <- 130
# c$g[grep("85",c$rcp)] <- 0
# c$b <- NA
# c$b[grep("26",c$rcp)] <- 255
# c$b[grep("45",c$rcp)] <- 255
# c$b[grep("60",c$rcp)] <- 45
# c$b[grep("85",c$rcp)] <- 0
#
#
# max.x <- max((UofC-x)/scalar)
#
# quartz(width = 10, height = 3, type = "pdf", "Figure5", file = "Figure5.pdf")
# par (fig=c(0,1,0,1), # Figure region in the device display region (x1,x2,y1,y2)
# 	omi=c(.6,.8,.3,.3), # global margins in inches (bottom, left, top, right)
# 	mai=c(0,.1,.1,0)) # subplot margins in inches (bottom, left, top, right)
# 	layout(matrix(c(1,2,3,4), 1, 4, byrow = TRUE), widths=c(3,3,3,3), heights=c(3))
#
# j <- c$rcp[1]
# y <- x[,grep(j,colnames(x))]
# dam <- data.frame(ECS=ECS, prob=probs, damages=(UofC-y[1,])/scalar) # Deterministic damages
# dam.det <- sample(dam$damages, size=10^6, replace=TRUE, prob=dam$prob)
# dam.det <- quantile(dam.det,seq(0,1,0.001))
# dam <- data.frame(ECS=rep(ECS,each=nrow(x)-1), prob=rep(probs,each=nrow(x)-1), damages=(UofC - as.vector(y[2:nrow(y),]))/scalar) # Stochastic damages
# dam.sto <- sample(dam$damages, size=10^6, replace=TRUE, prob=dam$prob)
# dam.sto <- quantile(dam.sto,seq(0,1,0.001))
#
# qqplot(dam.sto,dam.det, xlim=c(0,1.05*max.x), ylim=c(0,1.05*max.x), bty="n",main="",xaxs="i",yaxs="i",xlab="",ylab="",xaxt="n",yaxt="n",pch = 20,col=rgb((c$r[c$rcp==j]), (c$g[c$rcp==j]), (c$b[c$rcp==j]), max = 255, alpha = 255))
# lines(x=c(0,15*10^14), y=c(0,15*10^14))
# axis(side=1, at = seq(0,15*10^14,length.out=4), labels = c("",seq(5*10^14,15*10^14,length.out=3)/10^12))
# axis(side=2, at = seq(0,15*10^14,length.out=4), labels = seq(0,15*10^14,length.out=4)/10^12, las=2)
#
# # Risk premium
# y <- x1[,grep(j,colnames(x))]
# dam <- data.frame(ECS=ECS, prob=probs, damages=y[1,]) # Deterministic damages
# dam.det <- sum(dam$damages * (dam$prob/sum(dam$prob)))
# dam <- data.frame(ECS=rep(ECS,each=nrow(x)-1), prob=rep(probs,each=nrow(x)-1), damages=as.vector(y[2:nrow(y),])) # Stochastic damages
# dam.sto <- sum(dam$damages * (dam$prob/sum(dam$prob)))
# rp <- (dam.det-dam.sto)/scalar1
# text(x=0, y=0.8*15*10^14, labels="Risk premium",pos=4, cex=1.5)
# text(x=0, y=0.6*15*10^14, labels=paste(round(rp/10^12,0), "Tr", sep=""),pos=4, cex=1.5)
# text(x=0, y=0.4*15*10^14, labels=paste(round(100*rp/gdp,0), "%", sep=""),pos=4, cex=1.5)
#
# j <- c$rcp[2]
# y <- x[,grep(j,colnames(x))]
# dam <- data.frame(ECS=ECS, prob=probs, damages=(UofC-y[1,])/scalar) # Deterministic damages
# dam.det <- sample(dam$damages, size=10^6, replace=TRUE, prob=dam$prob)
# dam.det <- quantile(dam.det,seq(0,1,0.001))
# dam <- data.frame(ECS=rep(ECS,each=nrow(x)-1), prob=rep(probs,each=nrow(x)-1), damages=(UofC - as.vector(y[2:nrow(y),]))/scalar) # Stochastic damages
# dam.sto <- sample(dam$damages, size=10^6, replace=TRUE, prob=dam$prob)
# dam.sto <- quantile(dam.sto,seq(0,1,0.001))
#
# qqplot(dam.sto,dam.det, xlim=c(0,1.05*max.x), ylim=c(0,1.05*max.x), bty="n",main="",xaxs="i",yaxs="i",xlab="",ylab="",xaxt="n",yaxt="n",pch = 20,col=rgb((c$r[c$rcp==j]), (c$g[c$rcp==j]), (c$b[c$rcp==j]), max = 255, alpha = 255))
# lines(x=c(0,15*10^14), y=c(0,15*10^14))
# axis(side=1, at = seq(0,15*10^14,length.out=4), labels = c("",seq(5*10^14,15*10^14,length.out=3)/10^12))
# axis(side=2, at = seq(0,15*10^14,length.out=4), labels = rep("",4), las=2)
#
# # Risk premium
# y <- x1[,grep(j,colnames(x))]
# dam <- data.frame(ECS=ECS, prob=probs, damages=y[1,]) # Deterministic damages
# dam.det <- sum(dam$damages * (dam$prob/sum(dam$prob)))
# dam <- data.frame(ECS=rep(ECS,each=nrow(x)-1), prob=rep(probs,each=nrow(x)-1), damages=as.vector(y[2:nrow(y),])) # Stochastic damages
# dam.sto <- sum(dam$damages * (dam$prob/sum(dam$prob)))
# rp <- (dam.det-dam.sto)/scalar1
# text(x=0, y=0.6*15*10^14, labels=paste(round(rp/10^12,0), "Tr", sep=""),pos=4, cex=1.5)
# text(x=0, y=0.4*15*10^14, labels=paste(round(100*rp/gdp,0), "%", sep=""),pos=4, cex=1.5)
#
# j <- c$rcp[3]
# y <- x[,grep(j,colnames(x))]
# dam <- data.frame(ECS=ECS, prob=probs, damages=(UofC-y[1,])/scalar) # Deterministic damages
# dam.det <- sample(dam$damages, size=10^6, replace=TRUE, prob=dam$prob)
# dam.det <- quantile(dam.det,seq(0,1,0.001))
# dam <- data.frame(ECS=rep(ECS,each=nrow(x)-1), prob=rep(probs,each=nrow(x)-1), damages=(UofC - as.vector(y[2:nrow(y),]))/scalar) # Stochastic damages
# dam.sto <- sample(dam$damages, size=10^6, replace=TRUE, prob=dam$prob)
# dam.sto <- quantile(dam.sto,seq(0,1,0.001))
#
# qqplot(dam.sto,dam.det, xlim=c(0,1.05*max.x), ylim=c(0,1.05*max.x), bty="n",main="",xaxs="i",yaxs="i",xlab="",ylab="",xaxt="n",yaxt="n",pch = 20,col=rgb((c$r[c$rcp==j]), (c$g[c$rcp==j]), (c$b[c$rcp==j]), max = 255, alpha = 255))
# lines(x=c(0,15*10^14), y=c(0,15*10^14))
# axis(side=1, at = seq(0,15*10^14,length.out=4), labels = c("",seq(5*10^14,15*10^14,length.out=3)/10^12))
# axis(side=2, at = seq(0,15*10^14,length.out=4), labels = rep("",4), las=2)
#
# # Risk premium
# y <- x1[,grep(j,colnames(x))]
# dam <- data.frame(ECS=ECS, prob=probs, damages=y[1,]) # Deterministic damages
# dam.det <- sum(dam$damages * (dam$prob/sum(dam$prob)))
# dam <- data.frame(ECS=rep(ECS,each=nrow(x)-1), prob=rep(probs,each=nrow(x)-1), damages=as.vector(y[2:nrow(y),])) # Stochastic damages
# dam.sto <- sum(dam$damages * (dam$prob/sum(dam$prob)))
# rp <- (dam.det-dam.sto)/scalar1
# text(x=0, y=0.6*15*10^14, labels=paste(round(rp/10^12,0), "Tr", sep=""),pos=4, cex=1.5)
# text(x=0, y=0.4*15*10^14, labels=paste(round(100*rp/gdp,0), "%", sep=""),pos=4, cex=1.5)
#
# j <- c$rcp[4]
# y <- x[,grep(j,colnames(x))]
# dam <- data.frame(ECS=ECS, prob=probs, damages=(UofC-y[1,])/scalar) # Deterministic damages
# dam.det <- sample(dam$damages, size=10^6, replace=TRUE, prob=dam$prob)
# dam.det <- quantile(dam.det,seq(0,1,0.001))
# dam <- data.frame(ECS=rep(ECS,each=nrow(x)-1), prob=rep(probs,each=nrow(x)-1), damages=(UofC - as.vector(y[2:nrow(y),]))/scalar) # Stochastic damages
# dam.sto <- sample(dam$damages, size=10^6, replace=TRUE, prob=dam$prob)
# dam.sto <- quantile(dam.sto,seq(0,1,0.001))
#
# qqplot(dam.sto,dam.det, xlim=c(0,1.05*max.x), ylim=c(0,1.05*max.x), bty="n",main="",xaxs="i",yaxs="i",xlab="",ylab="",xaxt="n",yaxt="n",pch = 20,col=rgb((c$r[c$rcp==j]), (c$g[c$rcp==j]), (c$b[c$rcp==j]), max = 255, alpha = 255))
# lines(x=c(0,15*10^14), y=c(0,15*10^14))
# axis(side=1, at = seq(0,15*10^14,length.out=4), labels = c("",seq(5*10^14,15*10^14,length.out=3)/10^12))
# axis(side=2, at = seq(0,15*10^14,length.out=4), labels = rep("",4), las=2)
#
# # Risk premium
# y <- x1[,grep(j,colnames(x))]
# dam <- data.frame(ECS=ECS, prob=probs, damages=y[1,]) # Deterministic damages
# dam.det <- sum(dam$damages * (dam$prob/sum(dam$prob)))
# dam <- data.frame(ECS=rep(ECS,each=nrow(x)-1), prob=rep(probs,each=nrow(x)-1), damages=as.vector(y[2:nrow(y),])) # Stochastic damages
# dam.sto <- sum(dam$damages * (dam$prob/sum(dam$prob)))
# rp <- (dam.det-dam.sto)/scalar1
# text(x=0, y=0.6*15*10^14, labels=paste(round(rp/10^12,0), "Tr", sep=""),pos=4, cex=1.5)
# text(x=0, y=0.4*15*10^14, labels=paste(round(100*rp/gdp,0), "%", sep=""),pos=4, cex=1.5)
#
# mtext("Damages with", side=3, outer=T, line=0.2, cex=.9, at=0)
# mtext(expression(paste("deterministic ", Delta, "T", sep="")), side=3, outer=T, line=-1, cex=.9, at=0)
# mtext(expression(paste("Damages with stochastic ", Delta, "T (Trillion USD)", sep="")),side=1,outer=T, line=3.5)
#
# dev.off()

# Read in risk premia for aleatory uncertainty only (from "figure2.R") and for the risk interaction effect (the difference between the risk premia calculated above and the risk premia from "figure2.R")
d <- t(data.frame(aleatory=c(3,9,15,32), interaction=c(9, 25, 31, 46)))
colnames(d) <- c("RCP2.6","RCP4.5","RCP6.0","RCP8.5")

## IPCC AR5 colour scheme
c <- data.frame(rcp=c("26","45","60","85"))
c$r <- NA
c$r[grep("26",c$rcp)] <- 0
c$r[grep("45",c$rcp)] <- 121
c$r[grep("60",c$rcp)] <- 255
c$r[grep("85",c$rcp)] <- 255
c$g <- NA
c$g[grep("26",c$rcp)] <- 0
c$g[grep("45",c$rcp)] <- 188
c$g[grep("60",c$rcp)] <- 130
c$g[grep("85",c$rcp)] <- 0
c$b <- NA
c$b[grep("26",c$rcp)] <- 255
c$b[grep("45",c$rcp)] <- 255
c$b[grep("60",c$rcp)] <- 45
c$b[grep("85",c$rcp)] <- 0

quartz(width = 10, height = 3, type = "pdf", "Barchart", file = "Barchart.pdf")
par(mar=c(4,5,3,1))
barplot(d, horiz=TRUE, xlab="Risk premium (Trillion USD)", col=c("white"), yaxt="n", border="white", xlim=c(0,55))
for (j in 1:4) {
	barplot(c(rep(0,j-1),d[1,j]), horiz=TRUE, xlab="", xaxt="n", col=rgb((c$r[j]), (c$g[j]), (c$b[j]), max = 255, alpha = 255), border="white", add=TRUE)
	barplot(c(rep(0,j-1),d[2,j]), horiz=TRUE, xlab="", xaxt="n", col=rgb((c$r[j]), (c$g[j]), (c$b[j]), max = 255, alpha = 150), border="white", add=TRUE)
}
axis(side=2, labels=colnames(d), at=seq(0.7,4.3,length.out = 4), las=2, tick=FALSE)
text(x=-.3,y=seq(0.7,4.3,length.out = 4), paste(d[1,],"Tr",sep=""), col="white",pos=4)
text(x=d[1,]-.3,y=seq(0.7,4.3,length.out = 4), paste(d[2,]-d[1,],"Tr",sep=""), col="white",pos=4)
text(x=d[2,]-.3,y=seq(0.7,4.3,length.out = 4), paste(d[2,],"Tr (", round(d[2,]*100/80,0),"%)",sep=""), col="black",pos=4)

mtext(side=3, at=17, "Aleatory-only\nrisk premium", cex=.8)
mtext(side=3, at=39, "Aleatory x Epistemic\nrisk interaction effect", cex=.8)

dev.off()
